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ABSTRACT 

We present a new method based on a Bayesian hierarchical model to extract con- 
straints on cosmological parameters from SNIa data obtained with the SALT-II 
lightcurve fitter. We demonstrate with simulated data sets that our method delivers 
tighter statistical constraints on the cosmological parameters over 90% of the time, 
that it reduces statitical bias typically by a factor ^2 — 3 and that it has better 
coverage properties than the usual % 2 approach. As a further benefit, a full posterior 
probability distribution for the dispersion of the intrinsic magnitude of SNe is ob- 
tained. We apply this method to recent SNIa data, and by combining them with CMB 
and BAO data we obtain fi m = 0.28 ± 0.02, f2 A = 0.73 ± 0.01 (assuming w = -1) and 
O m = 0.28 ± 0.01, w — —0.90 ±0.05 (assuming flatness; statistical uncertainties only). 
We constrain the intrinsic dispersion of the B-band magnitude of the SNIa population, 
obtaining crjf* = 0.13 ± 0.01 [mag]. Applications to systematic uncertainties will be 
discussed in a forthcoming paper. 



Key words: Supernovae type la, Bayesian statistics, cosmological parameters, sys- 
tematic errors, intrinsic dispersion. 

1 INTRODUCTION 

Since the late 1990s when the Supernova Cosmology Project and the High-Z Supernova Search Team presented evidence that 

), observations of Type la supernovae 
(SNIa) has been seen as one of our most important tools for measuring the cosmic expansion as a function of time. Since 
a precise measurement of the evolution of the scale factor is likely the key to characterizing the dark energy or establishing 
that General Relativity must be modified on cosmological scales, the limited data that our universe affords us must be used 
to the greatest possible advantage. One important element of that task is the most careful possible statistical analysis of the 
data. Here we report how an improved statistical approach to SNIa data, with, in particular, more consistent treatment of 
uncertainties, leads to significant improvements in both the precision and accuracy of inferred cosmological parameters. 

The fundamental assumption underlying the past and proposed use of Type la supernovae to measure the expansion 
history is that they are "standardizable candles" . Type la supernovae occur when material accreting onto a white dwarf from 
a companion drives the mass of the white dwarf above the maximum that can be supported by electron degeneracy pressure, 
the Chandrasekhar limit of about 1.4 solar masses. This triggers the collapse of the star and the explosive onset of carbon 
fusion, in turn powering the supernova explosion. Because the collapse happens at a particular critical mass, all Type la 
supernovae are similar. Nevertheless, variability in several factors, including composition, rotation rate, and accretion rate, 
can lead to measurable differences in supernova observables as a function of time. Indeed, the intrinsic magnitude of the nearby 
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the expansion of the universe is accelerating (Riess et al.||1998 Perlmutter et al7||1999 
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Type la supernovae, the distances to which are known via independent means, exhibit a fairly large scatter. Fortunately, this 
scatter can be reduced by applying the so-called "Phillips corrections" - phenomenological correlations between the intrinsic 
magnitude of SNIa and their colour as well as between their intrinsic magnitudes and the time scale for the decline of their 
luminosity ( Phillips|1993 Phillips et al.|[l999[ |. Such corrections are derived from multi- wavelength observations of the SNIa 
lightcurves (i.e., their apparent brightness as a function of time). Fortuitously, they make SNIa into standardizable candles 
- in other words, having measured the colour and light curve of a SNIa, one can infer its intrinsic magnitude with relatively 
low scatter, typically in the range 0.1 — 0.2 mag. Observations of SNIa at a range of redshifts can then be used to measure 
the evolution of luminosity distance as a function of redshift, and thence infer the evolution of the scale factor, assuming that 
the intrinsic properties of SNIa do not themselves evolve (an assumption that has to be carefully checked). 

The SNIa sample which is used to measure distances in the Universe has grown massively thanks to a world-wide 

Wood-Vasey et al.fl2007 



observational effort (Astier et al. 2006 



2009a Freedman et al. 



2009 



Contreras et al. 



2010 



Amanullah et al. 2010 



Balland et al. 



2009} |Bailey et al. 



Kowalski et al. 2008 



2008 



Hicken et al. 



Kessler et al. 



20091. Presently, 



several hundred SNIa have been observed, a sample which is set to increase by an order of magnitude in the next 5 years 
or so. As observations have become more accurate and refined, discrepancies in their modeling have come into focus. Two 
main methods have emerged to perform the lightcurve fit and derive cosmological parameter constraints. The Multi-Colour 
Lightcurve Shape (MLCS) ( |Jha, Riess fc Kirshner|[2007 1 strategy is to simultaneously infer the Phillips corrections and the 
cosmological parameters of interest, applying a Bayesian prior to the parameter controlling extinction. The SALT and SALT- 
II ( Guy et al.,||2007 l methodology splits the process in two steps. First, Phillips corrections are derived from the lightcurve 
data; the cosmological parameters are then constrained in a separate inference step. As the supernova sample has grown and 
improved, the tension between the results of the two methods has increased. 

Despite the past and anticipated improvements in the supernova sample, the crucial inference step of deriving cosmological 
constraints from the SALT-II lightcurve fits has remained largely unchanged. For details of how the cosmology fitting is 
currently done, see for example (Astier et al. 2006 Kowalski et al. 2008 Amanullah et al. 2010 Conley et al. 20111. As 



currently used, it suffers from several shortcomings, such as not allowing for rigorous model checking, and not providing a 
rigorous framework for the evaluation of systematic uncertainties. The purpose of this paper is to introduce a statistically 
principled, rigorous, Bayesian hierarchical model for cosmological parameter fitting to SNIa data from SALT-II lightcurve 
fits. In particular the method addresses identified shortcomings of standard chi-squared approaches - notably, it properly 
accounts for the dependence of the errors in the distance moduli of the supernovae on fitted parameters. It also treats more 
carefully the uncertainty on the Phillips colour correction parameter, escaping the pontetial bias caused by the fact that the 
error is comparable to the width of the distribution of its value. 

We will show that our new method delivers considerably tighter statistical constraints on the parameters of interest, while 
giving a framework for the full propagation of systematic uncertainties to the final inferences. (This will be explored in an 
upcoming work.) We also apply our Bayesian hierarchical model to current SNIa data, and derive new cosmological constraints 
from them. We derive the intrinsic scatter in the SNIa absolute magnitude and obtain a statistically sound uncertainty on its 
value. 

This paper is organized as follows: in section [2] we review the standard method used to perform cosmological fits from 
SALT-II lightcurve results and we describe its limitations. We then present a new, fully Bayesian method, which we test on 
simulated data in section [3j where detailed comparisons of the performance of our new method with the standard approach 
are presented. We apply our new method to current SN data in section [4] and give our conclusions in section [5] 



2 COSMOLOGY FROM SALT-II LIGHTCURVE FITS 



2.1 Definition of the inference problem 



Several methods are available to fit SNe lightcurves, including the MLCS method, the Amis method, CMAGIC, (Wang et al 



2003 Conley et al., 2006) SALT, SALT-II and others. Recently, a sophisticated Bayesian hierarchical method to fit optical and 



infrared lightcurve data has been proposed by Mandel et al. (2009 20101. As mentioned above, MLCS fits the cosmological 



parameters at the same time as the parameters controlling the lightcurve fits. The SALT and SALT-II methods, on the 
contrary, first fit to the SNe lightcurves three parameters controlling the SN magnitude, the stretch and colour corrections. 
From those fits, the cosmological parameters are then fitted in a second, separate step. In this paper, we will consider the 
SALT-II method (although our discussion is equally applicable to SALT), and focus on the second step in the procedure, 
namely the extraction of cosmological parameters from the fitted lightcurves. We briefly summarize below the lightcurve 
fitting step, on which our method builds. 

The rest-frame flux at wavelength A and time t is fitted with the expression 

dF TCS t 



(IX 



-(t, A) = x [Mo (t, A) + x\M\ (t, A)] exp(c- CL{\)) . 



(1) 



where Mo, M\,CL are functions determined from a training process, while the fitted parameters are xo (which controls 
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the overall flux normalization), x\ (the stretch parameter) and c (the colour correction parameter). The .B-band apparent 
magnitude m* B is related to xq by the expression 



m B 



-2.5 log 



x / d\M Q (t = 0, A)T S (A) 



(2) 



where T B (A) is the transmission curve for the observer's .B-band, and t — is by convention the time of peak luminosity. After 



fitting the SNIa lightcurve data with SALT-II algorithm, e.g. Kessler et al. (2009a I report the best-fit values for m B ,X\,c 



the best-fit redshift z of each SNIa and a covariance matrix d for each SN, describing the covariance between m* B ,x\,c from 
the fit, of the form 



d = 



I 



2 



(77. 



2 



(3) 



' CI 



Let us denote the result of the SALT-II lightcurve fitting procedure for each SN as 

Di — {zi,m Bi ,xu,Ci,Ci}. (4) 

(where i runs through the n SNe in the sample, and measured quantities are denoted by a hat). We assume (as it is implicitly 
done in the literature) that the distribution of rh Bi ,xu,&i is a multi-normal Gaussian with covariance matrix d. 

The distance modulus fj,i for each SN (i.e., the difference between its apparent B-band magnitude and its absolute 
magnitude) is modeled as: 

Mi = rnht - Mi + a ■ xu ~ P ■ a (5) 

where Mi is the (unknown) .B-band absolute magnitude of the SN, while a, j3 are nuisance parameters controlling the stretch 
and colour correction (so-called "Phillips corrections"), respectively, which have to be determined from the data at the same 
time as the parameters of interest. The purpose of applying the Phillips corrections is to reduce the scatter in the distance 
modulus of the supernovae, so they can be used as almost standard candles. However, even after applying the corrections, some 
intrinsic dispersion in magnitude is expected to remain. Such intrinsic dispersion can have physical origin (e.g., host galaxies 
properties such as mass (jKelly et aL"||2010| |Sullivan et alj |2011| and star formation rate ( |Sullivan et al.|[2006[ ), host galaxy 



reddening ( Mandel et al.|20T0 i, possible SNe la evolution ( Gonzalez-Gaitan et al.|2011 1, etc) or be associated with undetected 
or underestimated systematic errors in the surveys. Below, we show how to include the intrinsic dispersion explicitly in the 
statistical model. 

Turning now to the theoretical predictions, the cosmological parameters we are interested in constraining are 

</f={Sl m ,Sl A or w,h} (6) 

where Sl m is the matter density (in units of the critical energy density), Qa is the dark energy density, w is the dark energy 
equation of state (taken to be redshift-independent, although this assumption can easily be generalized) and h is defined as 
Ho — 100/ikm/s/Mpc, where Ho is the value of the Hubble constant todajQ The curvature parameter fl K is related to the 
matter and dark energy densities by the constraint equation 

Q K = 1 - 0, n - fi A . (7) 

In the following, we shall consider either a Universe with non-zero curvature (with £7 K 7^ and an appropriate prior) but 
where the dark energy is assumed to be a cosmological constant, i.e. with w — —1 (the ACDM model), or a flat Universe 
where the effective equation of state parameter is allowed to depart from the cosmological constant value, i.e. Q K — 0, w 7^ — I 
(the wCDM model). 

In a Friedman-Robertson- Walker cosmology defined by the parameters ¥ ', the distance modulus to a SN at redshift Zi is 
given by 



m = \i[zi,f) = 5 log ^ p ' c " ' + 25. 

where Dl denotes the luminosity distance to the SN. This can be rewritten as 

fii = 77 + 5logd L (zi, Sl m , SIa, w), 

where 



7) 



-5 log 



100ft 



25 



(8) 



(9) 



(10) 



and c is the speed of light in km/s. We have defined the dimensionless luminosity distance (with Dl — c/HocLl, where c is 



1 The Hubble parameter h actually plays the role of a nuisance parameter, as it cannot be constrained by distance modulus measurements 
independently unless the absolute magnitude of the SNe is known, for the two quantities are perfectly degenerate. 
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the speed of light) 

d L (z, Sl m , 0\, w) = (1 jl^ sinn{y|aJ / dz' [(1 + z'ffl m + Q de (z) + (1 + z') 2 Q K ] ~ 1/2 } (11) 

Vitrei JO 

with the dark energy density parameter 

n dc (z) = Q A exp ^3 f* ^^r dx ) ■ ( 12 ) 

In the above equation, we have been completely general about the functional form of the dark energy equation of state, w(z). 
In the rest of this work, however, we will make the further assumption that w is constant with redshift, i.e., w(z) = w. We 
have defined the function sinn(a;) = x, sin(a;), sinh(x) for a flat Universe (fl K — 0), a closed Universe (fi K < 0) or an open 
Universe, respectively. 

The problem is now to infer, given data D in Eq. Q, the values (and uncertainties) of the cosmological parameters c ff, 
as well as the nuisance parameters {a, /3}, appearing in Eq. |5| and any further parameter describing the SNe population and 
its intrinsic scatter. Before building a full Bayesian hierarchical model to solve this problem, we briefly describe the usual 
approach and its shortcomings. 



2.2 Shortcomings of the usual \ method 



The usual analysis (e.g., Astier et al. (20061; Kowalski et al. (20081; Kessler et al. (2009al; Conley et al. ( 2011 1) defines a \ 2 
statistics as follows: 

0* ~ M? bs ) 2 



(13) 



where fii is given by Eq. |9| as a function of the cosmological parameters and the "observed" distance modulus fi° bs is obtained 
by replacing in Eq. (|5| the best-fit values for the colour and stretch correction and B-band magnitude from the SALT-II 
output (denoted by hats). Furthermore, the intrinsic magnitude for each SN, Mi, is replaced by a global parameter M, which 
represents the mean intrinsic magnitude of all SNe in the sample: 

/i° bs = m*Bi - M + a ■ xu - P ■ &i , (14) 

where the mean intrinsic magnitude M is unknown. The variance cr^ comprises several sources of uncertainty, which are 
added in quadrature: 

2 /fit\2,/2!\2,/int\2 
C/M = \ ff lti) + + \ a p ) > 

where af^ is the statistical uncertainty from the SALT-II lightcurve fit, 

2 



fit 



(15) 
(16) 



where ^ — (l,a,—/3) and Ci is the covariance matrix given in Eq. cr^j is the uncertainty on the SN redshift from 
spectroscopic measurements and peculiar velocities of and within the host galaxy; finally, <rjf is an unknown parameter 



describing the SN intrinsic dispersion. Further discussions of the unknown ctJJ 1 estimation problem, see ( Blondin,Mandel& 



Kirshner 



2011 



Kim 



2011 



Vishwakarma& Narlikar 



2011 1. As mentioned above, ideally crjf is a single quantity that encapsulates 



the remaining intrinsic dispersion in the SNIa sample, folding in all of the residual scatter due to physical effects not well 
captured by the Phillips corrections. However, observational uncertainties such as the estimation of photometric errors can 
lead to a variation of a 1 ^ sample by sample (for which there is a growing body of evidence). While we do not consider the 
latter scenario in this paper, it is important to keep in mind that describing the whole SN population with a single scatter 
parameter crjf is likely to be an oversimplification. 

Further error terms are added in quadrature to the total variance, describing uncertainties arising from dispersion due to 
lensing, Milky Way dust extinction, photometric zero-point calibration, etc. In this work, we do not deal with such systematic 
uncertainties, though they can be included in our method and we comment further on this below. 

The cosmological parameter fit proceeds by minimizing the x 2 in Eq. (131, simultaneously fitting the cosmological pa- 
rameters, q, P and the mean intrinsic SN magnitude M. The value of crjf* is adjusted to obtain x^/dof ~ 1 (usually on a 
sample-by-sample basis), often rejecting individual SNe with a residual pull larger than some arbitrarily chosen cut-off. It was 
recognized early that fitting the numerator and denominator of Eq. ( 13 1 iteratively leads to a large "bias" in the recovered 



value of P ( |Kowalski et al.||2008| |Astier et aL]|2006| |Wang et ai]|2006 



This has been traced back to the fact that the error 
on the colour correction parameter Cj is as large as or larger than the width of the distribution of values of Cj, especially 
for high-redshift SNe. This is a crucial observation, which constitutes the cornerstone of our Bayesian hierarchical model, 
as explained below. We demonstrate that an appropriate modeling of the distribution of values of d leads to an effective 



likelihood that replaces the \ °f Eq. ( 13 1. With this effective likelihood and appropriate Bayesian priors, all parameters can 
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be recovered without bias. If instead one adopts a properly normalized likelihood function, i.e., replacing the \ °f Eq. (13 1 
with 

(with the pre-factor Lq chosen so that the likelihood function integrates to unity in data space), marginalization over a, j3 
leads to catastrophic biases in the recovered values (up to ~ 6a in some numerical tests we performed). This is a strong hint 
that the naive form of the likelihood function above is incorrect for this problem. The effective likelihood we derive below 
solves this problem. 

The standard approach to cosmological parameters fitting outlined above adopted in most of the literature to date has 
several shortcomings, which can be summarized as follows: 

(i) The expression for the x 2 j Eq. ( |13| ), has no fundamental statistical justification, but is based on a heuristic derivation. 
The fundamental problem with Eq. ( |13[ ) is that some of the parameters being fitted (namely, a,/?) control both the location 
and the dispersion of the \ 2 expression, as they appear both in the numerator and the denominator, via the (er^*) 2 term. 
Therefore, the statistical problem is one of jointly estimating both the location and the variance. We show below how this 
can be tackled using a principled statistical approach. 

(ii) Adjusting a 1 ^ to obtain the desired goodness-of-fit is problematic, as it does not allow one to carry out any further 
goodness-of-fit test on the model itself, for obviously the variance has been adjusted to achieve a good fit by construction. 
This means that model checking is by construction not possible with this form of the likelihood function. 

(iii) It would be interesting to obtain not just a point estimate for ct™*, but an actual probability distribution for it, as 
well. This would allow consistency checks e.g. among different surveys, to verify whether the recovered intrinsic dispersions 
are mutually compatible (within errorbars). This is currently not possible given the standard y 2 method. A more easily 
generalizable approach is desirable, that would allow one to test the hypothesis of multiple SNe populations with different 
values of intrinsic dispersion, for example as a consequence of evolution with redshift, or correlated with host galaxy properties. 
Current practice is to split the full SN sample in subsamples (e.g., low- and high-redshift, or for different values of the colour 
correction) and check for the consistency of the recovered values from each of the subsamples. Our method allows for a more 
systematic approach to this kind of important model checking procedure. 

(iv) It is common in the literature to obtain inferences on the parameters of interest by minimizing (i.e., profiling) over 



nuisance parameters entering in Eq. ( 13 1. This is in general much more computationally costly than marginalization from e.g. 



MCMC samples ( |Feroz et alT]|2011[>. There are also examples where some nuisance parameters are marginalized over, while 



others are maximised (Astier et al. 20061, which is statistically inconsistent and should best be avoided. It should also be 
noted that maximisation and marginalization do not in general yield the same errors on the parameters of interest when the 
distribution is non-Gaussian. From a computational perspective, it would be advantageous to adopt a fully Bayesian method, 
which can be used in conjunction with fast and efficient MCMC and nested sampling techniques for the exploration of the 
parameter space. This would also allow one to adopt Bayesian model selection methods (which cannot currently be used with 
the standard \ 2 approach as they require the full marginalization of parameters to compute the Bayesian evidence). 



(v) The treatment of systematic errors is being given great attention in the recent literature (see e.g. Nordin et al. (20081) 



but the impact of various systematics on the final inference for the interesting parameters, has often been propagated in 



an approximate way (e.g., Kessler et al. (2009a I, Appendix F). As we are entering an epoch when SN cosmology is likely to be 
increasingly dominated by systematic errors, it would be desirable to have a consistent way to include sources of systematic 
uncertainties in the cosmological fit and to propagate the associated error consistently on the cosmological parameters. The 
inclusion in the analysis pipeline of systematic error parameters has been hampered so far by the fact that this increases 
the number of parameters being fitted above the limit of what current methods can handle. However, if a fully Bayesian 
expression for the likelihood function was available, one could then draw on the considerable power of Bayesian methods (such 
as MCMC, or nested sampling) which can efficiently handle larger parameter spaces. 

Motivated by the above problems and limitations of the current standard method, we now proceed to develop in the next 
section a fully Bayesian formalism from first principles, leading to the formulation of a new effective likelihood which will 
overcome, as will be shown below, the above problems. A more intuitive understanding of our procedure can be acquired from 
the simpler toy problem described in Appendix [B] 



2.3 The Bayesian hierarchical model 

We now turn to developing a Bayesian hierarchical model (BHM) for the SNe data from SALT-II lightcurve fits. The same 
general linear regression problem with unknown variance has been addressed by Kelly (20071, and applied in that paper 
to X-ray spectral slope fitting. The gist of our method is shown in the graphical network of Fig. [I] which displays the 
probabilistic and deterministic connection between variables. The fundamental idea is that we introduce a new layer of 
so-called "latent" variables - that is, quantities which describe the "true" value of the corresponding variables, and which 
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i (^y *(^) 




Figure 1. Graphical network showing the deterministic (dashed) and probabilistic (solid) connections between variables in our Bayesian 
hierarchical model (BHM). Variables of interest are in red, latent (unobserved) variables are in blue and observed data (denoted by hats) 
are in green. 



are obviously unobserved (represented in blue in Fig. [TJ . In particular, we associate with each SN a set of latent variables 
(zi, a,xu, Mi), which give the true value of the redshift, colour correction, stretch correction and intrinsic magnitude for that 
SN. For the stretch and colour correction, we have noisy estimates in the form of data Ci,xu (with an associated covariance 
matrix) resulting from the SALT-II fits (solid vertical lines). For the redshift Zi, we have noisy estimates given by the measured 
values Zi. As above, all observed quantities are denoted by a hat. 

Each of the latent variables (d, xh, Mi) is modeled as an uncorrelated random variable drawn from an underlying parent 
distribution, described by a Gaussian with two parameters, a population mean and a variance^] The latent intrinsic magnitudes 
Mi are thought of as realizations from a SN population with overall mean intrinsic magnitude Mo and intrinsic dispersion 
elf* (both to be estimated from the data): 



M t ~N{M Q ,{rf) 2 ). 



(19) 



This is represented in Fig. [I] by the solid arrow connecting the parameters Mo, (crjf*) 2 with the latent variables Mi. This 
represents statistically the physical statement that an intrinsic dispersion (u 1 ^) 2 is expected to remain in the distribution of 
SNe magnitudes, even after applying the Phillips corrections to the observed magnitudes. The simple Ansatz of Eq. ( |19| l can 
be replaced by a more refined modeling, which describes the existence of multiple SN populations, for example with intrinsic 
magnitude correlated with host galaxy properties, for which there is a growing body of evidence ( Sullivan et al.|200 6 Mandcl 
et al.|20l"0 Sullivan et al.pbll I. We shall investigate this possibility in a forthcoming work. 

The parent distribution of the true colour and stretch corrections is similarly represented by Gaussians, again parame- 
terised each by a mean (c^a;*) and a variance (R 2 ,R 2 ) as 

Cl ~ Af{c*, Rl), xu ~ Af{x*, Rl). (20) 

As above, the quantities c*, R 2 , x*, R 2 have to be estimated from the data. The choice of a Gaussian distribution for the 
latent variables c and x^ is justified by the fact that the observed distribution of c and x lt shown in Fig. [2] for the actual SNIa 
sample described in section |3.1| below, is fairly well described by a Gaussian. As shown in Fig. [2] there might be a hint for a 
heavier tail for positive values of c, but this does not fundamentally invalidate our Gaussian approximation. It would be easy 
to expand our method to consider other distributions, for example mixture models of Gaussians to describe a more complex 
population or a distribution with heavier tails, if non-Gaussianities in the observed distribution should make such modeling 



necessary. In this paper we consider the simple uni-modal Gaussians given by Eq. ( 20 \ 



While we are interested in the value (and distribution) of the SNe population parameters M , (crjf ) , the parameters 



2 In this paper we use the notation: 

x~Af(m,a 2 ) (18) 

to denote a random variable x being drawn from an underlying Gaussian distribution with mean m and variance a 2 . In vector notation, 
m is replaced by a vector m, while a 2 is replaced by the covariance matrix S and the distribution is understood to be the appropriate 
multidimensional Normal, see Appendix [A| 
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stretch, x x 



Figure 2. Histogram of observed stretch parameters xn and observed colour parameters 6j from the 288 SNIa from Kessler et al. 
(2009a), compared with a Gaussian fit (red curve). 



entering Eqs. (20 1 are not physically interesting. Therefore they will be marginalized out at the end. It is however crucial to 



introduce them explicitly in this step: it is precisely the lack of modeling in the distribution of that usually leads to the 
observed "biases" in the reconstruction of /3, see Appendix [B] for an explicit example with a toy model. For ease of notation, 



we re-write Eq. (191 and (201 in matrix notation as: 

<A/"(M ,E A ) 



M 



where 



(21) 
(22) 
(23) 

(24) 
(25) 

(26) 

Having introduced 3n latent (unobserved) variables (c,x lt M), where n is the number of SNe in the sample, the funda- 



■Af(c* ■ l„,diag(#e ■ 1„)) = p{c\c*,R a ) 
■N{Xi, ■ l n ,diag(i?^ • l n )) = p(x 1 \x i ,,R x ) 



M — (Mi, . . . , M n ) € 
M = Mo ■ l n 6 B", 

E A = diagf (o-JT*) 2 ■ ! 



mental strategy of our method is to link them to underlying population parameters via Eqs. (191 and (201, then to use the 



observed noisy estimates to infer constraints on the population parameters of interest (alongside the cosmological parameters), 
while marginalizing out the unobserved latent variables. 

The intrinsic magnitude Mi is related to the observed _B-band magnitude rh* B and the distance modulus /i by Eq. J5j, 
which can be rewritten in vector notation as: 

m* B = m + M- a?i + /3c. (27) 

Note that the above relation is exact, i.e. M, x_ x ,c are here the latent variables (not the observed quantities), while rn* B is the 
true value of the B-band magnitude (also unobserved). This is represented by the dotted (deterministic) arrows connecting 
the variables in Fig. [I] 

We seek to determine the posterior pdf for the parameters of interest O = {"W, a, /3, <xjf }, while marginalizing over 
the uknown population mean intrinsic magnitude, Mo. From Bayes theorem, the marginal posterior for is given by (see 
e.g. Trotta (20081; Hobson et al. (20101 for applications of Bayesian methods to cosmology): 

p(D \e,M )p(e,M ) 

p(D) 



dMn 



(28) 



P (e\D) = j dM oP (e,M \D) = 

where p(D) is the Bayesian evidence (a normalizing constant) and the prior p(Q, Mo) can be written as 

p(9, Mo) = p(V , a, p)p(M , <*) = p{Y, a, /3)p(M \a^)p(a^). (29) 

We take a uniform prior on the variables a, /3 (on a sufficiently large range so as to encompass the support of the likelihood), 
as well as a Gaussian prior for p(Mo|cJf ), since Mo is a location parameter of a Gaussian (conditional on cr™ ). Thus we take 



p{M 1 cr^ nt ) = JVmo (M m , <7m ) , 



(30) 



where the mean of the prior (M m = —19.3 mag) is taken to be a plausible value in the correct ballpark, and the variance 
(o"m = 2.0 mag) is sufficiently large so that the prior is very diffuse and non-informative (the precise choice of mean and 
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variance for this prior does not impact on our numerical results). Finally, the appropriate prior for cr™* is a Jeffreys' prior, 
i.e., uniform in logajf , as o™ is a scale parameter, see e.g. Box & Tiao (19921. Although the intrinsic magnitude Mo and 



the Hubble constant Ho are perfectly degenerate as far as SNIa data are concerned, we do not bundle them together in a 
single parameter but treat them separately with distinct priors, as we are interested in separating out the variability due to 
the distribution of the SNIa intrinsic magnitude. 

We now proceed to manipulate further the likelihood, p(D\Q, Mo) = p(c, x 1; m* B \0, Mo): 

p(c,%,ms|6,M ) = f dcd^ dMp(c,x 1 ,m* B \c,x 1 ,M,@,M )p(c,x 1 ,M\@,M ) (31) 
— J dc dx 1 dM p(c, x 1 , m* B \c, x lt M, €>) 

x j dR c dR* dc* da;* p(c|c*, R c )p(x 1 \x ir , R x )p{M\M , a i ^ t )p(R c )p(R x )p(c^)p(x i ,) (32) 

In the first line, we have introduced a set of 3n latent variables, {c, x 1; M_}, which describe the true value of the colour, 
stretch and intrinsic magnitude for each SNIa. Clearly, as those variables are unobserved, we need to marginalize over them. 
In the second line, we have replaced p(c, x t , MJO, Mo) by the distributions of the latent {c, x 1 , M} given by the probabilistic 



relationships of Eq. (211 and Eqs. (22 231, obtained by marginalizing out the population parameters {R c , Rx, c*, a:*}: 

p(c,x.i,M.\Q,Mo) = f dR c dR x dc* da;* p(cjc*, RJpfe^x*, R x )p(M\M , a^ t )p(R c )p(R x )p(c*)p(x ic ). (33) 



(we have also dropped Mo from the likelihood, as conditioning on Mo is irrelevant if the latent M are given) . If we further 
marginalize over Mo (as in Eq. ( |28[ ), including the prior on Mo), the expression for the effective likelihood, Eq. ( |32[ ), then 
becomes: 



p(c, x 1 , rh* B \Q) = j dc dx 1 dM p(c, x x , rn" B \c, a^, M, 6) 



x / dR c dR x dc* da;* dM p(c|c*, R c )p(x 1 \x ir , R x )p{M\M , a; lt )p(7? c )p(^)p(c*)p(a;*)p(Mo|o-jr 



(34) 

The term p(c, g 1; m* B \c, x lt M, Q) is the conditional probability of observing values {c^i^rnjg} if the latent (true) value 
of c, , M and of the other cosmological parameters were known. From Fig. [T] m* B is connected only deterministically to all 
other variables and parameters, via Eq. (271. Thus we can replace m* B = u + M — a • x x + • c and write 



p(c, m* B \c, x lt M, 6) = J^A/X^ii + Mi - a ■ xu + ft ■ a, Gi) 

i=l 

= |27rS c |-2 exp f -»[(* - Xa) T ^c\X - X Q ) 
where Hi = p,i(zi, O) and we have denned 

X = {Xi,. . .,X n } G R 3n , Xo = {Xo.i, . . . , X ,n} G M. 3 ", 
X z = {ci,xi :i , (Mi - ax^i + /3ci)} e R 3 , X , % = {cj, xi :i ,m Bi - /ij G K 3 , 
as well as the 3n x 3n block covariance matrbj^] 



E c = 



( dl 








^ 





<% 



















V o 








c n ) 



(35) 
(36) 

(37) 
(38) 

(39) 



Finally we explicitly include redshift uncertainties in our formalism. The observed apparent magnitude, rn* B , on the 



left-hand-side of Eq. (351, is the value at the observed redshift, z. However, /i in Eq. (351 should be evaluated at the true 
(unknown) redshift, z. As above, the redshift uncertainty is included by introducing the latent variables z and integrating 
over them: 



p(c,a: 1 ,M|c 1 a; 1 ,M, 9) = J dz p(c,x 1 ,M]c,x 1 ,M_,z,Q)p(z\z) 
where we model the redshift errors p(z\z) as Gaussians: 

z ~ J\f(z, E z ) 



(40) 
(41) 



3 Notice that we neglect correlations between different SNIa, which is reflected in the fact that Ec takes a block-diagonal form. It would 
be however very easy to add arbitrary cross-correlations to our formalism (e.g., coming from correlated systematic within survey, for 
example zero point calibration) by adding such non- block diagonal correlations to Eq. (|39|. 
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Parameter 


Symbol 


True Value 


Matter energy density parameter 




0.3 


Dark energy density parameter 


n A 


0.7 


Dark energy equation of state 


w 


-1 


Spatial curvature 


n K 


0.0 


Hubble expansion rate 


H [km/s/Mpc] 


72.0 


Mean absolute magnitude of SNc 


iWn fmael 

"-'U L oJ 


-19.3 


Intrinsic dispersion of SNe magnitude 


(j If 1 * [magi 


0.1 


Stretch correction 


a 


0.13 


Colour correction 




2.56 


Mean of distribution of x ± 




0.0 


Mean of distribution of c 




0.0 


s.d. of distribution of x^ 




1.0 


s.d. of distribution of c 


He 


0.1 


Observational noise on m* B 


"Si 


Depending on survey 


Observational noise on x x 


®x\ i 


Depending on survey 


Observational noise on c 




Depending on survey 


Correlation between x^ and c 




0.0 



Table 1. Input parameter values used in the generation of the simulated SNe SALT-II fits. 



with a n x n covariance matrix: 



E z = diag(cr zi , 



(42) 



It is now necessary to integrate out all latent variables and nuisance parameters from the expression for the likelihood, Eq. ( 34 1. 
This can be done analytically, as all necessary integral are Gaussian. The detailed steps are described in Appendix [C] In 
summary, the procedure consists of: 



(i) Marginalization over the intrinsic redshifts, see Eq. \C2\ . 

(ii) Marginalization over the latent variables {M,c, a^}, see Eq. ( C24 1 . 
(hi) Marginalization over the nuisance parameters {i*,c*}, see Eq. (C28l. 



This leads to the final expression for the effective likelihood of Eq. ( C32 1 : 

p(c,Xi,m B \Q) = I dlog.R c dlogiZ* |27rEcr2|27rEj»rs|27rE A |3|27rE r5|27rii'|! 



x exp 



hxfz^Xo - a t e a a - fcjirto + £eq h 



where the various vectors and matrices appearing in the above expression are defined in Appendix [C] This equation is the 
main result of our paper. The two remaining nuisance parameters R C ,R X cannot be integrated out analytically, so they need 
to be marginalized numerically. Hence, they are added to our parameters of interest and are sampled over numerically, and 
then marginalized out from the joint posterior. 



3 NUMERICAL TESTS OF THE METHOD 



3.1 Generation of simulated SNe data 

We now proceed to test the performance of our method on simulated data sets, and to compare it with the usual \ 2 approach. 



We take as starting point the recent compilation of 288 SNIa from Kessler et al. (2009a I, the data set on which we apply our 
Bayesian hierarchical model in the next section. Kessler et al. ( 2009a[ ) re-fitted the lightcurves from five different surveys: 



SDSS: 103 SNe ( [Kessler et al.|2009a[ ), 

ESSENCE: 56 SNe ( jMiknaitis fe Pignata|2007| |Wood-Vasey et al.|2007) , 
SNLS: 62 SNe ( |Astier et al.|2006l ), 

Nearby Sample: 33 SNe ( |Jha, Riess fc Kirshner|2007| > and 

HST: 34 SNe ( |Garnavich et al.|1998| |Knop et al.|2003)|Riess fc Strolger||2004| [2007] ), 



using both the SALT-II method and the MLCS fitter. In the following, we are exclusively employing the results of their 
SALT-II fits and use those as the observed data set for the purposes of our current work, as described in the previous section. 
More refined procedures could be adopted, for example by simulating lightcurves from scratch, using e.g. the publicly available 
package SNANA ( Kessler et al.|200 9b). In this paper we chose a simpler approach, consisting of simulating SALT-II fit results 
in such a way to broadly match the distributions and characteristics of the real data set used in Kessler et al. (2009a I. 
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Figure 3. An example realization of our simulated data sets (coloured according to survey), superimposed on real data (black). Colour 
code for simulated data survey: nearby sample (cyan), ESSENCE (green), SDSS (red), SNLS (blue) and HST (magenta). 



The numerical values of the parameters used in the simulation are shown in Table [T] We adopt a vanilla, flat ACDM 
cosmological model as fiducial cosmology. The Phillips correction parameters are chosen to match the best-fit values reported 



in Kessler et al. ( 2009a I , while the distributional properties of the colour and stretch correction match the observed distribution 
of their total SN sample. For each survey, we generate a number of SNe matching the observed sample, and we model their 
redshift distribution as a Gaussian, with mean and variance estimated from the observed distribution within each survey. The 
observational error of m* B ,c, is again drawn from a Gaussian distribution whose mean and variance have been matched to 
the observed ones for each survey. Finally, the pseudo-data (i.e., the simulated SALT-II fits results) are generated by drawing 
from the appropriate distributions centered around the latent variables. For simplicity, we have set to the off-diagonal 



elements in the correlation matrix ( 39 1 in our simulated data, and neglected redshift errors. None of these assumptions have 



a significant impact on our results. In summary, our procedure for each survey is as follows: 

(i) Draw a value for the latent redshift Zi from a normal distribution with mean and variance matching the observed ones. 
As we neglect redshift errors in the pseudo-data for simplicity (since the uncertainty in z is subdominant in the overall error 
budget), we set Zi = Zi. 

(ii) Compute fit using the fiducial values for the cosmological parameters "af and the above z% from Eq. Q. 

(iii) Draw the latent parameters xu,Ci,Mi from their respective distributions (in particular, including an intrinsic scatter 
0/f = 0-1 ma g in the generation of Mi). 

(iv) Compute m Bi using xu,c%, Mi and the Phillips relation Eq. JHJ. 

(v) Draw the value of the standard deviations <T xl i,<T 0i ,a mi , from the appropriate normal distributions for each survey 
type. A small, ^-dependent stochastic linear addition is also made to a Xl i, a Ci , <r mi , to mimic the observed correlation between 
redshift and error. 

(vi) Draw the SALT-II fit results from xu ~ N(xu, cr xi i), &i ~ M{ci, a Ci ) and rh Bi ~ Af(m Bi , o mi ). 

As shown in Fig. [3j the simulated data from our procedure have broadly similar distributions to the to real ones. The 
two notable exceptions are the overall vertical shift observed in the distance modulus plot, and the fact that our synthetic 
data cannot reproduce the few outliers with large values of the variances (bottom panels). The former is a consequence of the 
different intrinsic magnitude used in our simulated data (as the true one is unknown). However, the intrinsic magnitude is 
marginalized over at the end, so this difference has no impact on our inferences. The absence of outliers is a consequence of 
the fact that our simulation is a pure phenomenological description of the data, hence it cannot encapsulate such fine details. 
While in principle we could perform outlier detection with dedicated Bayesian procedures, we do not pursue this issue further 
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Parameter ACDM wCDM 

Uniform: W(0.0, 1.0) Uniform: W(0.0,1.0) 
fl K Uniform: W(-1.0, 1.0) Fixed: 

w Fixed: -1 Uniform: U{— 4, 0) 

Hp [km/s/Mpc] Af(72,8 2 ) Af(72,8 2 ) 

Common priors 
o-jf' [mag] Uniform on logo-jf*: W(-3. 0,0.0) 

M [mag] Uniform: «(-20.3, -18.3) 

a Uniform: 1/(0.0, 1.0) 

/8 Uniform: U(0.0, 4.0) 

ii c Uniform on log i? c : «(-5.0, 2.0) 

ij x Uniform on log R c : U(-5.0, 2.0) 



Table 2. Priors on our model's parameters used when evaluating the posterior distribution. Ranges for the uniform priors have been 
chosen so as to generously bracket plausible values of the corresponding quantities. 



in this paper. We stress once more that the purpose of our simulations is not to obtain realistic SNIa data. Instead, they 
should only provide us with useful mock data sets coming from a known model so that we can test our procedure. More 
sophisticated tests based on more realistically generated data (e.g., from SNANA) are left for future work. 



3.2 Numerical sampling 



After analytical marginalization of the latent variables, we are left with the following 8 parameters entering the effective 



likelihood of Eq. (431 



{Vim, fi« or w, H , cr^, a, f3, R c , R x } . (43) 

As mentioned above, in keeping with the literature we only consider either flat Universes with a possible w 7^ — 1 (the ACDM 
model), or curved Universes with a cosmological constant (w — —1, the wCDM model). Of course it is possible to relax those 
assumptions and consider more complicated cosmologies with a larger number of free parameters if one so wishes (notably 
including evolution in the dark energy equation of state). 

Of the parameters listed in Eq. (431, the quantities R C ,R X are of no interest and will be marginalized over. As for the 



remaining parameters, we are interested in building their marginal 1 and 2-dimensional posterior distributions. This is done by 
plugging the likelihood ( 43 1 into the posterior of Eq. ( 28 1 , with priors on the parameters chosen according to Table [2] We use a 
Gaussian prior on the Hubble parameter Ho = 72 ±8 km/s/Mpc from local determinations of the Hubble const ant ( [Freedman 
et al.|200T k However, as Ho is degenerate with the intrinsic population absolute magnitude Mo (which is marginalized over 
at the end), replacing this Gaussian prior with a less informative prior I/o [km/s/Mpc] ~ W(20, 100) has no influence on our 
results. 

Numerical sampling of the posterior is carried out via a nested sampling algorithm (Skilling 2004, 2006. Feroz & Hobson 



2008 Feroz et al. 20091. Although the original motivation for nested sampling was to compute the Bayesian evidence, the 



recent development of the MultiNest algorithm ( Feroz & Hobson 2008 ; Feroz et al.|2009 1 has delivered an extremely powerful 
and versatile algorithm that has been demonstrated to be able to deal with extremely complex likelihood surfaces in hun- 
dreds of dimensions exhibiting multiple peaks. As samples from the posterior are generated as a by-product of the evidence 
computation, nested sampling can also be used to obtain parameter constraints in the same run as computing the Bayesian 
evidence. In this paper we adopt the publicly available MultiNest algorithm (Feroz & Hobson 2008) to obtain samples from 
the posterior distribution of Eq. ( |28[ ). We use 4000 live points and a tolerance parameter 0.1, resulting in about 8 x 10 5 
likelihood evaluations^ 

We also wish to compare the performance of our BHM with the usually adopted x 2 minimization procedure. To this end, 



we fit the pseudo-data using the x expression of Eq. ( 13 1. In order to mimic what is done in the literature as closely as possible, 
we first fix a value of crjf*. Then, we simultaneously minimize the \ 2 w.r.t. the fit parameters ■& — {Q m , 0« or w, Ho, Mo, a, /3}, 
as described below. We then evaluate the x 2 /dof from the resulting best fit point, and we adjust cjf* to obtain x 2 /dof = 1. 
We then repeat the above minimization over ■& for this new value of crJJ 1 , and iterate the procedure. Once we have obtained 
the global best fit point, we derive 1- and 2-dimensional confidence intervals on the parameters by profiling (i.e., maximising 
over the other parameters) over the likelihood 



L(i?) = exp 



(44) 



with x given by Eq. (13 1. According to Wilks' theorem, approximate confidence intervals are obtained from the profile 



4 A Fortran code implementing our method is available from the authors upon request. 
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Figure 4. Reconstruction of cosmological parameters from a simulated data set encompassing 288 SNIa, with characteristics matching 
presently available surveys (including realization noise). Blue regions contain 95% and 68% of the posterior probability (other parameters 
marginalized over) from our BHM method, the red contours delimit 95% and 68% confidence intervals from the standard \ 2 method 
(other parameters maximised). The yellow star indicates the true value of the parameters. The left panel assumes w = — 1 while the 
right panel assumes Q K = 0. Notice how out method produced considerably less biassed constraints on the parameters. 



likelihood as the regions where the \ 2 increases by A\ 2 from its minimum value, where A\ 2 can be computed from the 
chi-square distribution with the number of degree of freedoms corresponding to the number of parameters of interest and is 
given in standard look-up tables. 

Obtaining reliable estimates of the profile likelihood using Bayesian algorithms (such as MultiNest) is a considerably harder 
numerical task than mapping out the Bayesian posterior. However, it has been shown that MultiNest can be successfully used 
for this task even in highly challenging situations ( Feroz et al.|20lT |, provided the number of live points and tolerance value 
used are adjusted appropriately. For our % 2 scan, we adopt 10 4 live points and a tolerance of 0.1. We have found that those 
values give accurate estimates of the profile likelihood more than 2a into the tails of the distribution for an 8 dimensional 
Gaussian toy model (whose dimensionality matches the case of interest here). With these MultiNest settings, we gather 
1.5 x 10 5 samples, from which the profile likelihood is derived. 

Our implementation of the \ 2 method is designed to match the main features of the fitting procedure usually adopted in 
the literature (namely, maximisation of the likelihood rather than marginalization of the posterior, and iterative determination 
of the intrinsic dispersion), although we do not expect that it exactly reproduces the results obtained by any specific imple- 
mentation. Its main purpose is to offer a useful benchmark against which to compare the performance of our new Bayesian 
methodology. 



3.3 Parameter reconstruction 

We compare the cosmological parameters reconstructed from the standard % 2 method and our Bayesian approach in Fig. [| for 
a typical data realization. The left-hand-side panel shows constraints in the fi m — CIa plane for the ACDM model, both from 
our Bayesian method (filled regions, marginalized 68% and 95% posterior) and from the standard x 2 method (red contours, 
68% and 95% confidence regions from the profile likelihood). In the right-hand-side panel, constraints are shown in the w — Q m 
plane instead (this fit assumes a flat Universe). In a typical reconstruction, our Bayesian method produced considerably tighter 
constraints on the cosmological parameters of interest than the usual \ 2 approach. Our constraints are also less biassed w.r.t. 
the true value of the parameters, an important advantage that we further characterize below. 

Our BHM further produces marginalized posterior distributions for all the other parameters of the fit, including the 
stretch and colour corrections and the intrinsic dispersion of the SNe. The ID marginal posteriors for those quantities are 
shown in Fig. [5] The recovered posterior means lie within 1<t of the true values. Notice that we do not expect the posterior 
mean to match exactly the true value, because of realization noise in the pseudo-data. However, as shown below, our method 
delivers less biassed estimates of the parameters, and a reduced mean squared error compared with the standard \ 2 approach. 
The stretch correction a is determined with 8% accuracy, while the colour correction parameter /3 is constrained with an 
accuracy better than 3%. A new feature of our method is that it produces a posterior distribution for the SN population 
intrinsic dispersion, crjf* (right-hand-side panel of Fig [Bj|. This allows to determined the intrinsic dispersion of the SNIa 
population to typically about 10% accuracy. 
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Figure 5. Marginalised posterior for the stretch correction a, colour correction parameter and logarithm of the intrinsic dispersion of 
SNe, logcrJJ 1 *, from a simulated data set from our Bayesian method. The vertical, dashed line gives the true value for each quantity. 



3.4 Comparison of long-term performance of the two methods 

Encouraged by the good performance of our BHM on single instances of simulated data, we now wish to compare the long-term 
performance of the two methods for a series of simulated data realizations. We are interested in comparing the average ability 
of both methods to recover parameter values that are as much as possible unbiased with respect to their true values, as well 
as to establish the coverage properties of the credible and confidence intervals. 

Coverage is defined as the probability that an interval contains (covers) the true value of a parameter, in a long series 
of repeated measurements. The defining property of a e.g. 95% frequentist confidence interval is that it should cover the 
true value 95% of the time; thus, it is reasonable to check if the intervals have the properties they claim. Coverage is a 
frequentist concept: intervals based on Bayesian techniques are meant to contain a given amount of posterior probability for 
a single measurement (with no reference to repeated measurements) and are referred to as credible intervals to emphasize the 
difference in concept. While Bayesian techniques are not designed with coverage as a goal, it is still meaningful to investigate 
their coverage properties. To our knowledge, the coverage properties of even the standard x 2 method (which, being a frequentist 
method would ideally be expected to exhibit exact coverage) have never been investigated in the SN literature. 

We generate 100 realizations of the pseudo-data from the fiducial model of Table [l] as described in section |3.1| and we 
analyze them using our BHM method and the standard \ 2 approach, using the same priors as above, given in Table [2] For 
each parameter of interest 0, we begin by considering the relative size of the posterior 68.3% range from out method, af HM , 
compared with the 68.3% confidence interval from the \ 2 method, ag , which is summarized by the quantity Sg which shows 
the percentage change in errorbar size with respect to the errorbar derived using the \ 2 method 



BHM 



x 100. 



(45) 



A value Se < 1 means that our method delivers tighter errorbars on the parameter 9. A histogram of this quantity for the 
variables of interest is shown in Fig. [6] from which we conclude that our method gives smaller errobars on Q m , SI a and w in 
almost all cases. The flip side of this result is that the uncertainty on a, /3 is larger from our method than from the x 2 approach 
in essentially all data realization. This is a consequence of the large number of latent variables our method marginalizes over. 

Tight errorbars are good, but not if they come at the expense of a biased reconstruction. To investigate this aspect, we 
build the following test statistics from each reconstruction: 



l/#truc — 1| — l^ys/^truc _ 1|, 



(46) 



where #bhm is the posterior mean recovered using our BHM method, 8 x2 is the best-fit value for the parameter recovered 
using the standard x 2 approach and 8 t ruc is the true value for that parameter. The meaning of Te is the following: for a given 
data realization, if the reconstructed posterior mean from our BHM is closer to the true parameter value than the best-fit 
X 2 , then Te < 0, which means that our method is less biassed than \ 2 ■ A histogram of the distribution of Te across the 100 
realizations, shown in Fig. [7J can be used to compare the two methods: a negative average in the histogram means that the 
BHM outperforms the usual x 2 . For all of the parameters considered, our method is clearly much less biassed than the \ 2 
method, outperforming x 2 about 2/3 of the time. Furthermore, the reconstruction of the intrinsic dispersion is better with 
our Bayesian method almost 3 times out of 4. We emphasize once more that our methodology also provides an estimate of 
the uncertainty in the intrinsic dispersion, not just a best-fit value as the x 2 approach. 

We can further quantify the improvement in the statistical reconstruction by looking at the bias and mean squared error 
(MSE) for each parameter, defined as 



bias = (9 — <9t ruc } 
MSE = bias 2 + Var, 



(47) 
(48) 
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Table 3. Comparison of the bias and mean squared error for our Bayesian method and the usual x 2 approach. The columns labelled 
"Improvement" give the factor by which our Bayesian method reduces the bias and the MSE w.r.t. the \ 2 approach. 
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Figure 6. Histograms of the quantity defined in Eq. (45) , comparing of the errorbars on each parameter from our method and from the 
standard \ 2 approach for 100 realization, for the ACDM model (left) and the wCDM model (right). A change in errorbar size of —10% 
indicates BHM errorbars are 10% smaller than x 2 errorbars. A change in errorbar size of +10% indicates BHM errorbars are 10% larger 
than x 2 errorbars. Our BHM method generally delivers smaller errors on the cosmological parameters (top row) but larger errors on the 
Phillips correction parameters (bottom row). 



respectively, where the expectation is taken by averaging over the observed values in our 100 simulated trials, 9 = #bhm 
(8 — 9^ 2 ) for the BHM (for the \ 2 approach) and Var is the observed parameter variance. The bias is the expectation value 
of the difference between estimator and true value, while the MSE measures the average of the squares of the errors, i.e., the 
amount by which the estimator differs from the true value for each parameter. Obviously, a smaller bias and a smaller MSE 
imply a better performance of the method. The results for the two methods are summarized in Table [3j which shows how 
our method reduces the bias by a factor ~ 2 — 3 for most parameters, while reducing the MSE by a factor of ~ 2. The only 
notable exception is the bias of the EOS parameter w, which is larger in our method than in the x 2 approach. 

Finally, in Fig. [8] we plot the coverage of each method for 68% and 95% intervals. Errobars give an estimate of the 
uncertainty of the coverage result, by giving the binomial sampling error from the finite number of realizations considered, 
evaluated from the binomial variance as Np(l — p), where N = 100 is the number of trials and p is the observed fractional 
coverage. Both method slightly undercover, i.e. the credible region and confidence intervals are too short, although the lack 
of coverage is not dramatic: e.g., the typical coverage of the la (2a) intervals from our method is ~ 60% (90%). Our method 
shows slightly better coverage properties than the \ 2 method, while producing considerably tighter and less biassed constraints 
(as demonstrated above). This further proves that the tighter intervals recovered by our method do not suffer from bias w.r.t 
the true values. 



4 COSMOLOGICAL CONSTRAINTS FROM CURRENT SNIA DATA 



We now apply our BHM to fitting real SN data. We use the SALT-II fits result for 288 SNIa from |Kessler et al.| ( |2009a| ), which 
have been derived from 5 different surveys. Our method only includes statistical errors according to the procedure described in 
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Figure 7. Histograms of the test statistics defined in Eq. | |46[ ), comparing the long-term performance of the two methods for the 
parameters of interest in the ACDM model (left) and the wCDM model (right). A predominantly negative value of the test statistics 
means that our method gives a parameter reconstruction that is closer to the true value than the usual x 2 > i.e., less biassed. For the 
cosmological parameters (top row), our method outperforms \ 2 about 2 times out of 3. 
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Figure 8. Coverage of our method (blue) and standard \ 2 (red) for 68% (solid) and 95% (dashed) intervals, from 100 realizations 
of pseudo-data for the ACDM model (left) and the wCDM model (right). While both methods show significant undercoverage for all 
parameters, our method has a comparable coverage to the standard x 2 , except for w. Coverage values for the intrinsic dispersion (T^ nt 
are not available from the \ 2 method, as it does not produce an error estimate for this quantity. 



section [2~3| coming from redshift uncertainties (arising from spectroscopic errors and peculiar velocities), intrinsic dispersion 
(which is determined from the data) and full error propagation of the SALT-II fit results. Systematic uncertainties play an 
important role in SNIa cosmology fitting, and (although not included in this study) can also be treated in our formalism in a 
fully consistent way. We comment on this aspect further below, though we leave a complete exploration of systematics with 
our BHM to a future, dedicated work (Ma rch et al.poTT l. 

We show in Fig. [9] the constraints on the cosmological parameters VL m — £7a (left panel, assuming w = — 1) and w — Q. m 
(right panel, assuming fl K = 0) obtained with our method. All other parameters have been marginalized over. In order to 
be consistent with the literature, we have taken a non-informative prior on Hq, uniform in the range [20, 100] km/s/Mpc. 
The figure also compares our results with the statistical contours from Kessler et al. (2009a I, obtained using the \ 2 method. 

we combine 



10 



(Notice that we compare with the contours including only statistical uncertainties for consistency.) In Fig. 
our SNIa constraints with Cosmic Microwave Background (CMB) data from WMAP 5-yrs measurements (Komatsu et al 



2009 



2005 



and Baryonic Acoustic Oscillations (BAO) constraints from the Sloan Digital Sky Survey LRG sample ( Eisenstein et al 



^ - r. 



, using the same method as Kessler et al. (2009a I. The combined SNIa, CMB and BAO statistical constraints result in 
: 0.28 ± 0.02, Sl A = 0.73 ±0.01 (for the ACDM model) and fi m = 0.28 ±0.01, w = -0.90 ±0.05 (68.3% credible intervals) 
for the wCDM model. Although the statistical uncertainties are comparable to the results by Kessler et al. (2009a) from the 
same sample, our posterior mean values present shifts of up to ~ 2a compared to the results obtained using the standard 
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Figure 9. Constraints on the cosmological parameters Q m , f^A (left panel, assuming to = —1) and w, Q m (right panel, assuming Q K = 0) 
from our Bayesian method (light/dark blue regions, 68% and 95% marginalized posterior), compared with the statistical errors from the 
usual x 2 approach (yellow/red regions, same significance level; from Kessler at al. (2009a)). The yellow star gives the posterior mean 
from our analysis. 



X approach. This is a fairly significant shift, which can be attributed to our improved statistical method, which exhibits a 
reduced bias w.r.t. the % 2 approach. 

Tig. [Tl] shows the Id marginalized posterior distributions for the Phillips correction parameters and for the intrinsic 
dispersion. All parameters are well constrained by the posterior, and we find a — 0.12 ±0.02, f3 = 2.7 ±0.1 and a value of the 
intrinsic dispersion (for the whole sample) a 1 ^ = 0.13 ± 0.01 mag. 



Kessler et al. (2009a I find values for the intrinsic dispersion ranging from 0.08 (for SDSS-II) to 0.23 (for the HST sample) 



but their \ 2 method does not allow them to derive an error on those determinations. With our method, it would be easy to 
derive constraints on the intrinsic dispersion of each survey - all one needs to do is to replace Eq. ( |19[ ) with a corresponding 
expression for each survey. This introduces one pair of population parameters (Mo, crjf') for each survey. In the same way, one 
could study whether the intrinsic dispersion evolves with redshift. We leave a detailed study of these issues to a future work. 
The value of a found in Kessler et al. (2009a I is in the range 0.10 — 0.12, depending of the details of the assumptions 

0.015. These results are comparable with our own. As for the colour 
vary in the range 2.46 



made, with a typical statistical uncertainty of order 



correction parameter /3, constraints from Kessler et al 
of order 0.1 



(2009a 



2.66, with a statistical uncertainty 



0.2. This stronger dependence on the details of the analysis seems to point to a larger impact of systematic 
uncertainties for /3, which is confirmed by evidence of evolution with redshift of the value of (3 (Kessler et al. (2009a), Fig. 39). 
Our method can be employed to carry out a rigorous assessment of the evolution with redshift of colour corrections. A possible 
strategy would be to replace /3 with a vector of parameters /3i,/?2, • • • , with each element describing the colour correction in 
a different redshift bin. The analysis proceeds then as above, and it produces posterior distributions for the components of 
P, which allows to check the hypothesis of evolution. Finally, in such an analysis the marginalized constraints on all other 
parameters (including the cosmological parameters of interest) would automatically include the full uncertainty propagation 
from the colour correction evolution, without the need for further ad hoc inflation of the errorbars. These kind of tests will 
be pursued in a forthcoming publication (March et al.|2011). 



5 CONCLUSIONS 

We have presented a statistically principled approach for the rigorous analysis of SALT-II SNIa lightcurve fits, based on a 
Bayesian hierarchical model. The main novelty of our method is that it produces an effective likelihood that propagates uncer- 
tainties in a fully consistent way. We have introduced an explicit statistical modeling of the intrinsic magnitude distribution of 
the SNIa population, which for the first time allows one to derive a full posterior distribution of the SNIa intrinsic dispersion. 

We have tested our method using simulated data sets and found that it compares favourably with the standard x 2 
approach, both on individual data realizations and in the long term performance. Statistical constraints on cosmological 
parameters are significantly improved, while in a series of 100 simulated data sets our method outperforms the \ 2 approach 
at least 2 times out of 3 for the parameters of interest. We have also demonstrated that our method is less biassed and has 
better coverage properties than the usual approach. 

We applied our methodology to a sample of 288 SNIa from multiple surveys. We find that the flat ACDM model is 
still in good agreement with the data, even under our improved analysis. However, the posterior mean for the cosmological 
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Figure 10. Combined constraints on the cosmological parameters Q m ,Q^ (left panel, assuming ui = —1) and w,Q m (right panel, 
assuming Q K = 0) from SNIa, CMB and BAO data. Red contours give 68% and 95% regions from CMB alone, green contours from BAO 
alone, blue contours from SNIa alone from our Bayesian method. The filled regions delimit 68% and 95% combined constraints, with the 
yellow star indicating the posterior mean. 




Figure 11. Marginalised posterior for the stretch correction a, colour correction ft parameter and logarithm of the intrinsic dispersion 
of SNe, log (r^ nt from current SNIa data. 



parameters exhibit up to 2a shifts w.r.t. results obtained with the conventional \ 2 approach. This is a consequence of our 
improved statistical analysis, which benefits from a reduced bias in estimating the parameters. 

While in this paper we have only discussed statistical constraints, our method offers a new, fully consistent way of including 
systematic uncertainties in the fit. As our method is fully Bayesian, it can be used in conjunction with fast and efficient Bayesian 
sampling algorithms, such as MCMC and nested sampling. This will allow to enlarge the number of parameters controlling 
systematic effects that can be included in the analysis, thus taking SNIa cosmological parameter fitting to a new level of 
statistical sophistication. The power of our method as applied to systematic errors analysis will be presented in a forthcoming, 
dedicated work. 

At a time when SNIa constraints are entering a new level of precision, and with a tenfold increase in the sample size 
expected over the next few years, we believe it is timely to upgrade the cosmological data analysis pipeline in order to extract 
the most information from present and upcoming SNIa data. This paper represents a first step in this direction. 
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APPENDIX A: NOTATION 



For reference, we collect here a few useful formulas relating to Gaussian integrals. Use the notation x ~ Nx_([i, E) to denote 
that the random variable x is drawn from a Normal distribution of mean fi and inverse covariance matrix E, given by 



p(x) 



1 



exp 



(x — /i) T E 1 (x — fi) 



In performing Gaussian integrals, it is also convenient to recall that: 

A/kQ^, Ei) • A4(m 2 , E 2 ) = fo ■ A4(/£ /; £/) 



where 



/o=AT Hi (/x 2 ,(E 1 +E 2 )) 

E / = (sr 1 + s 2 - 1 )- 1 . 



(Al) 

(A2) 

(A3) 
(A4) 
(A5) 



Finally, the integral of the multidimensional Gaussian of Eq. ( Al I over all space is unity. 



APPENDIX B: TOY MODEL FOR GENERAL LINEAR FITTING 



We can get an intuitive understanding of the central ingredient in our Bayesian hierarchical method by considering a simpler 
toy model, which highlights the salient features of the problem. 

As mentioned in section [2^2] several authors have observed in the past that the spread of fit values for the colour correction, 
c, is of the same order as the size of the statistical uncertainty on the values themselves. This leads to a bias in the best-fit 



value of P obtained from minimizing the x i n Eq. ( |13[ ). The quantity /3 gives the slope of the linear relationship between c 
and fx, see Eq. (271. In the usual x' 2 approach, the latent (true) c are replaced by the observed value, c, as in Eq. (141. If 



the statistical uncertainty in the independent variable, c, is as large as the spread of its values, the best-fit slope obtained 
from minimizing the x 2 mav be biased, as the large uncertainty in the actual location of the independent variable leads to 
confusion as to the value of the slope. The (Bayesian) solution to this problem is to determine the spread of the independent 
variable directly from the data, and to marginalize over it with an appropriate prior. This gives the most general solution to 
the problem of linear fitting with errors in both the independent and dependent variable, as shown by Gull ( 19891). Recent 
literature in cosmology and astronomy ( D'Agostini|2 005 ; Hog g et al.|2010 1 addresses linear fitting, but not this general case, 
which has been treated before in Kelly (20071 (see also Andreon (20061; Andreon & Hurn (20101 for related examples). In 



this short appendix, we will analyse the toy model of fitting a linear function, and compare the performance of a BHM and 
the x 2 approach. 



Bl Bayesian linear fitting in the presence of large x and y uncertainties 



In this subsection, we give a short review of the results of Gull ( 1989 1. The simplest toy model which illustrates the methodology 



we adopt in our paper is shown by the graphical network of Fig. |B1| A linear relationship is assumed between the latent 
variables Xi and yi, described by a slope a and intercept b (which are collectively denoted as a parameter vector 6): 



yi = axi + b. 

The observed values for the dependent and independent variables are denoted by hats (xi,yi, i = 1 
obtained from the latent values under the assumption of Gaussian noise (with known variances a 2 , a 2 . 



(Bl) 

, N), and they are 



Xi ~ Af(xi,al) and yi ~ M{yi, a 2 ) 



(B2) 
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Figure Bl. Bayesian network showing the deterministic and probabilistic connections in the toy linear model. Solid lines indicate 
probabilistic connections, dashed lines represent deterministic connections. Parameters to be constrained are in red, latent variables in 
blue and data in green (denoted by hats). 9 denotes the parameters a, b, i.e. the intercept and slope of the linear relation of Eq. JbTI. 



This probabilistic relationship is depicted in Fig. |B1| by the solid arrows connecting the latent variables to the observed 
quantities. Assuming that errors are uncorrelated, the joint likelihood is given bjj^] 



p{x,y\x,y,a x ,o y ,6) = (4tt a x a y )~ 



exp 



+ 



£i(& 



Vi 



The problem can be made more symmetric by defining rescaled versions of the data: 



X 



X — Xo 
Rx 



and Y 



y-yo 

Ry 



(B3) 



(B4) 



where the variables xo, yo describe the mean value of x, y, while R x , R y describe their spread. These new variables are related 
to the old ones (a, b) by 

b — yo — axo and a = R y /R x . (B5) 
Neglecting the normalization constant, the joint posterior for x, y, xo,yo, Rx, R y can be written as: 

p{x,y,x ,yo,Rx,R y \x,y,a x ,a y ) oc p(x, y\x, x Q , y , Rx, R y )p(x\x , yo, Rx, R y )p(x , y , Rx, Ry), (B6) 



where the first term on the r.h.s. is the likelihood of Eq. (B3|. The key step is to recognize that the appropriate conditional 
distribution for the latent x is 



p(x\x ,yo,R x ,R y ) = (4tt 2 RI) 



-1/2 



exp 



2 



x ) 



m 



(B7) 



This describe a prior over x centered around the hyperparameter xo and with standard deviation Rx. Crucially, both xo and 
R x are unknown, and are explicitly determined from the data in the joint posterior, before being marginalized out at the 



end. Finally, for the prior p(xq, yo, R x , R y ) appearing in Eq. (B6l we adopt a uniform prior on xo,yo (as those are location 
variables) and a prior uniform in log R x , log R y (those being scale variables, as apparent from (|B7l). 



From here, one can eliminate y from Eq. (B6l using Eq. (Bl I, then trade 6 for Ry using Eq. (B5|, to obtain 

2 



p(x,xo,yo,a,logR x ,logRy\x,y,tj x ,a y ) oc (87r a x a y R 



R 2N-iV/2 



exp 



+ 



yo + ax ) Y,ii x i - x o) 



Rl 
(B8) 

From this expression, the latent x can be marginalized out analytically, as well as the nuisance parameters xo,yo, by using 
appropriate completions of the square in the Gaussian. After some algebra, one finally obtains 



p(a,logR\x,y,a x ,a y ) (x (a 2 a x Rl + a x al + alR x ) 2 e xp 



1 V xx {a 2 R 2 x + al) - 2V xy aR 2 x + V vy {R 2 x + a 2 ) 



where: 



V xx — ^ ^ { Xi 



I a 2 alR x + ct§ct| + a^R x 

v xy = ^2(&i-x)(yi-y), V 2 y = ^ - y) 2 , R=(R x R y ) 1/2 . 



(B9) 
(BIO) 



For ease of notation, quantities without subscript i denote in the following iV-dimensional vectors, e.g. x = {x\, . . . ,x^}. 
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Simulated data x 2 Constraints 




L °9io a L °9io a 

Figure B2. Numerical comparison between the x 2 approach and the Bayesian method for linear fitting. Upper left panel: data set of 
N = 300 observations with errors both in the x and y directions, given by the error bar. Upper right panel: reconstruction of the slope a 
and intercept b using a x 2 likelihood (red cross is the true value, green circle the maximum likelihood value). Lower left panel: Bayesian 
posterior (Eq. | |B9[ |) in the log a, log it! plane, with green circle showing posterior mean. In both panels, contours enclose la, 2a and 
3<r regions. Lower right panel: marginalized Bayesian posterior (blue) and profiled x 2 likelihood (black, lying exactly on top of the blue 
curve), with the dashed line showing the true value. The two methods give essentially identical results in this case. 



In the above, x — £\ Xi/N, and similarly for y. The marginal posterior for the slope a is obtained by numerical marginalization 
of log 7? from the above expression. 



B2 Comparison with the \ 2 approach 

If instead of the statistically principled solution found above, one writes down a simple \ 2 expression for the likelihood 
including error propagation from the linear relationship ( |B1[ ) one would obtain (D'Agostini |2005"| ): 

_b) 2 

2 



L ^^\-\T, \7a^' (BU) 
from where the intercept 6 is eliminated by maximising over it (profiling). We now compare the reconstruction of the slope 



parameter a from Eq. (Bill (with b eliminated via profiling) with the result obtained using the Bayesian expression, (B9l 



marginalized numerically over log R (with a uniform prior on log R) with the help of simulated data. 

Figs. |B2| and |B3| show in the upper left panel the simulated data points (N = 300), with the error bars giving the size of the 
standard deviation in the x and y directions for each datum (i.e., the value of a x , a y ). The contour plots depict joint posterior 



regions for (log a, log R) from the Bayesian expression of Eq. (B9l, and confidence regions from the \ 2 likelihood (Bll I in the 
a, b plane. The lower right panel shows the Id marginalized posterior distribution for log a (blue) and the profile likelihood 
from the x 2 approach (black). Fig. B2 shows that the two results are essentially identical when the size of the statistical 
error is smaller than the spread of values of the data points (in this particular example, by a factor of 2). However, when the 
statistical uncertainty in the x direction is as large as or larger than the spread of the points (Fig. B3 1 , the x 2 approach gives 



a biassed result for the slope parameter a. The Bayesian posterior, on the contrary, is closer to the true value, although it 
(correctly) shows a larger uncertainty. 



APPENDIX C: DERIVATION OF THE EFFECTIVE LIKELIHOOD 
CI Integration over intrinsic redshifts 

In order to perform the multi-dimensional integral over z, we Taylor expand fi around £ (as justified by the fact that redshift 
errors are typically small: the error from 300 km/s peculiar velocity is a Zi = 0.0012, while the error from spectroscopic redshifts 
from SNe themselves is a Zi — 0.005, see Kessler et al. ( 2009a l): 

'D L {zj) 



5 log^ 



Mpc 



+ 25«/z(%) + 5(log 10 e) 



d Zj D L (z 3 ) 



(CI) 



With this approximation we can now carry out the multi-dimensional integral of Eq. (40 1, obtaining 
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Figure B3. As in Fig. |B3| but now with a larger statistical uncertainty in the data compared to their spread. The Bayesian marginal 
posterior for the slope (blue, lower right panel) peaks much closer to the true value than the x 2 expression (black). In the lower left 
panel, the 3cr contour from the Bayesian method lies outside the range of the plot. 



p(m* B \c, x 1; M, Q) = \2nS m \ * 

1 



x exp 



(ZBs — (m + M - c)) T S m 1 (m| 3 — (p, + M — a ■ x 1 + ■ c)) 



where from now on, /x = fj,(z) and 

/ = diag(/i,...,/„) 
D' L (zi) 



ft = 51og 10 (e) 
_ 51og 10 (e) 



D L (z 
D L (z,) 



1+Zi 



+ &i) x cosn{ v / r a^[ f dz {{i + z') 3 n m + n dc (z) + (i + z ) 2 n K ] 1/2 } 

Jo 



(C2) 

(C3) 
(C4) 

(C5) 



(C6) 
(C7) 



D L { Zi ) 

x((i + z'fn m + n dc (z) + (i + z) 2 fi K n 1/2 )] 

Strictly speaking, one should integrate over redshift in the range ^ z% < oo, not — oo < zi < oo, which would result in the 
appearance of Gamma functions in the final result. However, as long as ^ <C 1 (as is the case here), this approximation is 
expected to be excellent. 



C2 Integration over latent jc, x 1; M \ 



From Eq. ( 32 l and using the expression in Eq. ( 35 1 , we wish to integrate out the latent variables 

Y = {Yi,...,Y n } GE 3n , 
Yi = {ci,x 1:il Mi} G R 3 , 



We therefore recast expression ( 35 I as 



KMi,m^|c,aL 1 ,M,e) = |E c r 5 exp( - -[(AY - Xo) 1 T,^ (AY - X ) 



where we have defined the block-diagonal matrix 



A = diag(T,T,...,T) e: 



with 





1 





" 




Ci 


T = 





1 







Xi 




_ P 


—a 


1 




Mi 



(08) 
(C9) 
(CIO) 

(en) 

(C12) 
(C13) 
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Improved constraints from SNIa 23 

The prior terms appearing in Eq. (34l, namely p(c\c t , R c )p(x 1 \x+, R x )p(M\Mq, cr.'f*), may be written as: 
p(c\c^R c )p(x 1 \x^R m )p(M\M ,cr^) = \2nE P \~i exp (~[(Y - Y.) T Hp\Y - Y,)\ 



where 



S 1 = diag(7? c 2 ,R X 2 , (cr™*) 2 



Ep 1 =diag(5 _1 ) S ,_1 ,...,S ,_1 ) G 



y, = J ■ b e 



' 1 


" 


1 








1 


1 





1 





. o 


1 _ 


c* 




x* 


e K 


M 





(C14) 

(C15) 
(C16) 

(C17) 



b = 

Now the integral over AY — Ac da^ dM in Eq. ( 34 1 can be performed, giving 
dY p(c, , | c, , M, Q)p(c| e* , i? c 1 x+ , R x )p(M\ M , cr™* ) = 

= |27rE c r'|27rEp|"2|27rEA|^ exp 

where 

v,-l aT^ 



(C18) 



I 1 = A T T,r 1 A + Ep 1 S R 3 " x3 ", 



E^Fo = A^E^Xo + Ep 1 !;, 

Fo = Ea(A t E5 1 Zo + Ep 1 F)E A (A + Ep 1 ^), 



A = A T T,a 1 X G R 3nxl . 



(C19) 

(C20) 
(C21) 
(C22) 
(C23) 



Substituting Eq. ( |C19[ ) back into Eq. (34| gives: 

wIb|9) = J d-R c di?* dc* dx* |27rSc| - ' |27rEp|"2 |27rE 4 |2 

x exp (-^[xZz^Xo - Y T T,- A l Y + YjZ^Y* 
x p(R c )p(R x )p(c*)p(x*)p(M \al?). 



(C24) 



C3 Integration over population variables {c*, i*, Mo} 



The priors on the population variables b = {c*, a;*, Mo} in Eq. ( C24 1 can be written as: 

P(k) = p(c ir )p(x*)p(M \a™ t ) = |2^S |-i exp (-\{b - 6 m f Eq \b - bjj 



where 



and 



— 



o 






1/a 



2 



b m = 






M„ 



(C25) 



(C26) 



(C27) 
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Thus Eq. (C24l can be written as: 



p{c,Xi,m B \Q) = J dR c dR x db |27rEcr 5 |27rEpr2|27rE. 4 | 5 |27rE r5p(7? c )p(^) 

x exp (-~[xZEc 1 X - (E A (A + Ep 1 J ■ 6)) T E2 1 (E A (A + Ep 1 J • 6)) + b T J T E P L Jb + (b- b m ) T ^\b - b„ 
= y dR c dRx \2nTtc\~^2irZp\-i\2nVA\~^2irV \~h(Rc)p(Rx) 

x exp (-^[XfE^-Xo - A T E A A - fejif-^o + 6^0 X &jJ j Ab exp f-i[(6 - k ) T K-\b - fc )M 



(C28) 



where 



K' 1 = — J t E p 1 EaEp 1 J + J T 'Sp 1 J + Eq 1 £ K 3x3 , (C29) 
K^ko = J^E^EaA + E,7 l b m e R 3x \ (C30) 
ko = ^(^Ep'EaA + Eo 1 ^). (C31) 
We can now carry out the Gaussian integral over 6 in Eq. ( |C24[ ), obtaining our final expression for the effective likelihood, 

P{c,x 1 ,m%\e)= [ dlog.R c dlogi4 \2irHc\"^ \2nY, P \~^ |2ttEa| ' ^TrEops \2nK\^ 



x exp ^-i^Eo'Xo - A t EaA - k^K^ko + gF^bj) , (C32) 
where we have chosen an improper Jeffreys' prior for the scale variables R c , R x : 

p(Rc) oc R" 1 => p(R c )dR c oc dlog7? c , (C33) 

and analogously for R x . These two remaining nuisance parameters cannot be integrated out analytically, so they need to be 
marginalized numerically. Hence, R c , R x are added to our parameters of interest and are sampled over numerically, and then 
marginalized out from the joint posterior. 
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